A new estimation method for COVID-19 time-varying reproduction number using active cases

We propose a new method to estimate the time-varying effective (or instantaneous) reproduction number of the novel coronavirus disease (COVID-19). The method is based on a discrete-time stochastic augmented compartmental model that describes the virus transmission. A two-stage estimation method, which combines the Extended Kalman Filter (EKF) to estimate the reported state variables (active and removed cases) and a low pass filter based on a rational transfer function to remove short term fluctuations of the reported cases, is used with case uncertainties that are assumed to follow a Gaussian distribution. Our method does not require information regarding serial intervals, which makes the estimation procedure simpler without reducing the quality of the estimate. We show that the proposed method is comparable to common approaches, e.g., age-structured and new cases based sequential Bayesian models. We also apply it to COVID-19 cases in the Scandinavian countries: Denmark, Sweden, and Norway, where the positive rates were below 5\% recommended by WHO.


Introduction
The coronavirus disease 2019 (COVID-19), a disease outbreak of atypical pneumonia that originated from Wuhan, China [1], has caused globally at least 250 million confirmed cases, including an estimated 5 million deaths in approximately 221 countries and territories by November 2021. The World Health Organization (WHO) declared the COVID-19 crisis a pandemic on 11 March 2020.
In modelling the disease's transmission as well as to inform and evaluate control policies, it is particularly important to estimate its reproduction number.
Early estimates for COVID-19 basic reproduction number R 0 , that denotes the transmission potential of infectious disease when introduced to a completely susceptible population, ranged from 1.4 to 6.49 [2]. The effective (or instantaneous) reproduction number R t , on the other hand, reflects the extent of transmission in the presence of population immunity or intervention. Thus, the estimation of R t is important for evaluating public measure success. However, estimation of R t is sensitive to the model structure and parameter assumptions [3]. As a case in point, due to incorporation of more individual case information and travel data, the estimate for R 0 in Wuhan was revised upward from 2.2-2.7 to 5.7 [4]. On the other hand, data inavailability or poor quality often hinders the use of certain estimation methods, such as serial interval data that are usually needed to estimate R t (e.g., Fraser [5], Wallinga and Teunis [6], Cauchemez et al. [7], White and Pagano [8]).
In the course of calculating the exact value of R t , especially when the data has not yet reached its peak, precise assumptions and data estimates are needed. Nishiura et al. [9] discussed a likelihood-based approach to estimate R t from early epidemic growth data, while Cazelles at al. [10] used stochastic models for the disease dynamics coupled with particle Markov chain Monte Carlo algorithm. Using the compartmental Susceptible-Infectious-Recovered (SIR) model, Bettencourt and Ribeiro [11] use the incidence data to estimate R 0 and R t . In this paper, based on the Susceptible-Infectious-Recovered-Dead (SIRD) model as a reference, we develop a novel approach to estimate R t of COVID-19. It uses information on the number of infected or active (I), recovered (R), and death (D) cases, which are readily available for all affected countries, so that they can be accessed rather easily. This method does not require information regarding serial intervals, which makes the estimation procedure simpler without reducing the quality of the estimate. We assume mass population testing is sufficiently enough, such that the positive rate is below 5% recommended by WHO. This is to ensure data quality of the number of infection is acceptable since asymptomatic carrier transmission is often underestimated [12].
The reproduction number is estimated from reported cases under uncertainties using a two-stage estimation method based on the Extended Kalman Filter (EKF) and a low-pass filter. The method not only considers the nominal number of reported cases, but also its daily pattern. To show our method's practical ability, we apply it to COVID-19 cases in the Scandinavian countries: Denmark, Sweden, and Norway, and compare the results with two commonly used Bayesian methods due to Bettencourt and Ribeiro [11] and Cori et al. [5,13]. We show that the results are indeed comparable. Remark that a similar approach, developed independently, can be found in [14]. The difference is in computational technique to estimate the reproduction number. In this paper, we estimate the reproduction number using EKF, while in [14] it was estimated using Kalman smoother.

A Discrete-Time Stochastic Augmented Compartmental Model
Our estimation method is based on the compartmental SIRD model that can be written as the following first-order nonlinear differential equations: where S, I, R, and D denote the number of susceptible cases, the number of active cases, the number of recovered cases, and the number of deceased cases, respectively. N is the total number of population, β is the average number of contacts per person per time, while γ and κ are the recovery and death rate.
Remark that the value of β is time-varying due to intervention, i.e., β=β(t). To use the model, we require information on the average infectious time T i and the Case Fatality Rate (CFR), so that For COVID-19, we take T i = 9 as the infectious period on average lasts for 9 days (7-11 days with 95% CI) [15], while the CFR is assumed around 1%. The time-varying effective reproduction number is then given by: The approximation is under the assumption that government intervention is taken at an early stage so that the susceptible is relatively the same over time as the total population. This is the case especially for emerging diseases. We modify the SIRD model by augmenting the following two equations into the system:Ė The former equation takes into account the daily number of new reported cases E, while the latter one says that the effective reproduction number R t is assumed to be a piece-wise constant function with jump every one day time interval.
Discretizing the model using the forward Euler method, we obtain the fol-lowing discrete-time augmented SIRD model: Our method computes a new estimate of R t based on new reported cases. Since their frequency is low (could be once a day), the reported data can be interpolated using, e.g., a modified Akima cubic Hermite interpolation, such that it fits with the time step ∆t. In our simulation, the time step ∆t is chosen as 0.01, i.e., 100 time discretization within one day interval. The confidence interval of our estimated R t is determined by computing the reproduction number for different values of the infectious period T i within a certain interval.
To simplify the presentation, we define the augmented state vector and as such, the discrete-time augmented SIRD model (8)- (13) can be written as follows where f is the nonlinear term written in the right hand side of (8)-(13) and w is introduced as an uncertainty to model the inaccuracies due to simplification in the modelling. The uncertainty is assumed to be a zero mean Gaussian white noise with known covariance Q F . This is to simplify the calculation since the actual epidemic data usually follow Gamma distribution. In practice, Q F can be considered as a tuning parameter for the EKF. Thus, the transmission model becomes a discrete-time stochastic augmented SIRD model. Reported cases, such as the number of active cases and the cumulative numbers of recovered and death, can be incorporated into the model using the following output vector Here, v denotes uncertainties due to false testing results. We also assume the uncertainty to be a zero mean Gaussian white noise with known covariance R F .
As well as Q F , R F can also be considered as a tuning parameter. Following the available data that include I, R, D, and E, the data/measurement matrix C is taken to be A Two-Stage Filtering Method A two-stage filtering method is used to estimate the daily reproduction number R t . The method consists of the EKF and a low-pass filter. In the first stage of estimation, the EKF is used to estimate the state variables and the value of R t under uncertainties in the number of reported cases. Afterwards, the low pass filter is used to remove short term fluctuations of the reported cases that can be caused by delays in the reporting. For example, suddenly in Denmark there were 893 recovered patients reported on 1 April 2020, in contrast to the previous days from 16 February 2020 onwards when there was no recovery reported at all. Such an accumulated delay can cause a falsely decreasing value of R t .
The EKF is an extension of Kalman filter for nonlinear systems. The Kalman filter itself is based on a recursive Bayesian estimation and is an optimal linear filter. The idea of EKF is to linearize the non-linearity around its estimate. Due to that linearization, the optimality and stability of the EKF cannot be guaranteed. However, if the non-linearity is not severe, the EKF can give a reasonably good estimate.
Let us denotex(k) as an estimated vector state from the EKF. Applying first-order Taylor series expansion to f atx(k), we obtain where J f (x(k)) is the Jacobian matrix of f , given by: where The EKF consists of two steps: predict and update. The discrete-time stochastic augmented SIRD model is used to predict the next state and covariance and update them after obtaining new data/measurement. The EKF can be considered as one of the simplest dynamic Bayesian networks. While the EKF calculates estimates of the true values of states recursively over time using incoming measurements and a mathematical process model, recursive Bayesian estimation calculates estimates of an unknown probability density function recursively over time using incoming measurements and a mathematical process model [16]. Letx(n|m) denotes the estimate of x at time n given observations up to and including at time m ≤ n. The Kalman filter algorithm is given as follows [17] Predictx Updateỹ (k + 1) = y(k + 1) − Cx(k + 1|k) (28) x(k + 1|k + 1) =x(k + 1|k) + K(k + 1)ỹ(k + 1) (30) Here P (k|k) denotes a posteriori estimate covariance matrix. In the second stage, a low pass filter based on a rational transfer function is used to remove short term fluctuation at time step k, and is given bŷ where y n is a window length along the data. In our case, we choose y n = 3 ∆t .
To evaluate the quality of the estimate, we calculate a Relative Root Mean Square Error (RRMSE) between the estimated and reported cases. The RRMSE is defined as where N d is the number of observed days and X ∈ {I, D, R, E}.
Case study: Scandinavian countries  In this section, we apply our method to study viral transmission of COVID-19 in Denmark, Sweden, and Norway. All datasets and MATLAB code are available on GitHub (https://github.com/agusisma/covid19). As of January 2021, the three Scandinavian countries have higher cumulative testing rate compared to other parts of the world, with Denmark held a record with 260 tests per 1000 population. During that time the daily test-positivity rate is below 1% for Denmark and Norway, while Sweden is at 2.9% [18]. These numbers are good indications about the testing capacity in the Scandinavian countries and may describe the dynamics of the transmission better with respect to asymptomatic cases. The countries also have a different approach in their public measures in responding to COVID-19, e.g., Sweden did not implement a strict lockdown, unlike its Nordic neighbouring countries.
We plot the observed incidence of COVID-19 in Denmark, Sweden, and Norway in Fig. 1. We also plot in the same figure estimated numbers computed using our method, where good agreement is obtained. For all estimation, the process and observation covariance matrices are considered as tuning parameters and are chosen as Q F = diag(10 10 10 10 5 0.2) and R F = diag(100 10 10 5 1), respectively. These parameters are obtained from trial and error and are chosen such that the RRMSE between the estimated and reported data are sufficiently small. In our case study, the RRMSE are shown in Table 1. Here, we can observe the method provides relatively small estimation errors for all countries. Norway has the largest error, which can be attributed to the lack of daily update of the active and recovered cases. In our simulation, we use the same tuning parameters. The error can be reduced by using different value of Q F and R F .
In applying our method, we also compare it with two commonly used methods to estimate transmission parameters, namely the sequential Bayesian method of Bettencourt and Ribeiro [11] that provides an approximation of the basic reproduction number, and the instantaneous method by Fraser [5] that is implemented with a Bayesian analysis [19,13]. The former method exploits the new reported incidence, while the latter one uses the distribution of the serial interval.  First, we compare our method with Bettencourt and Ribeiro [11], that allows sequential estimation of the basic reproduction number at the initial stage when the growth is still exponential. While the two methods are based on the SIR model, Bettencourt and Ribeiro [11] use new incidence data and the result is filtered using a five-day moving average filter. In Fig. 3, we plot the comparison and summarise the basic reproduction numbers that are taken to be the maximum of the curves in Table 2. It is interesting to note how the methods give rather similar estimations. This indicates that our method gives comparable results to those of [11].
Finally, we plot the time-varying effective reproduction number in Fig. 4. Here, we compare our results with those using Cori et al. [13]. The method of [13] utilises the disease serial interval, which we approximate using a shifted Gamma distribution [13] with mean 4.7 and standard deviation 2.9 [20]. The prior belief for the value of R t is taken to be Gamma function with mean and standard deviation 5. We do not average out the data of daily new cases, but instead take the likelihood estimation of a new case at one day to depend also on the estimation of the previous three days.
In Fig. 4 we obtain that the two methods give the plot of R t with the same trend, indicating that our method is also comparable with [13]. There is a delay of about four days in the trend, especially with the time when the reproduction number curve crossed the horizontal axis. The delay is caused by the peaks of new daily cases and active ones that also differ by about the same days.
A different trend especially at later times between the methods appears in   Fig. 4(c) for Norway. The curve from our method is quite smooth, while it is rather fluctuating in that using Cori et al. [13]. The discrepancy is caused by the active and recovered cases that apparently were not updated regularly, in contrast to the new positive cases needed by the method of [13]. The unreported recovery cases were all released at once on May 22, 2020, see Fig. 1(c).

Conclusion and Future Work
Many mathematical models and estimation methods have been developed to estimate several types of reproduction numbers during epidemic outbreaks. Here, we provide a novel method exploiting reported active, recovered and death cases using the SIR model as an underlying approach. This new method offers several advantages compared to existing methods: (i) from modeling point of view, the resulting R t value can follow the dynamics of the model suggested, so it is possible to develop it further if the model chosen has a higher complexity, (ii) the estimation method can still be expanded in terms of statistical view, and (iii) the method does not need information about serial intervals. In the case that the data provided in time series do not change much or instead have drastic changes, such as accumulating at a certain time, the resulting R t value will show the same spikes and serrations. As a result, the latest information from data dynamics can be more elaborated.
By applying the method to COVID-19 cases in the Scandinavian countries and comparing the results to commonly used methods due to [11] and [13], we showed that our model is comparable, which expectedly will allow for fast assessment of the reproduction number in new outbreaks. Using the method to forecast and critically assess incidence data in countries with high underreporting, such as Indonesia, is addressed for future work.

Data and Code Availability
Data and codes can be found here: https://github.com/agusisma/covid19